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INTRODUCTION 


The objective of the present research is to develop a general mathemati- 
cal model and solution methodologies for analyzing structural response of 
thin, metallic shell-type structures under large transient, cyclic or static 
thermomechanical loads. Among the system responses, which are associated with 
these load conditions, are thermal buckling, creep buckling and racheting. 
Thus, geometric as well as material -type nonlinearities (of high order) can be 
anticipated and must be considered in the development of the mathematical 
model. Furthermore, this must also be accommodated in the solution procedures. 

SUMMARY OF PROGRESS 


The progress to date has been elaborated upon in an interim scientific 
report submitted to the sponsor during the summer of 1986, and in a series of 
semiannual progress reports. The most recent of these is dated April, 1988. 

A complete true ab-initio rate theory of kinematics and kinetics for 
continuum and curved thin structures, without any restriction on the magnitude 
of the strains or the deformation, was formulated. The time dependence and 
large strain behavior are incorporated through the introduction of the time 
rates of the metric and curvature in two coordinate systems; a fixed (spatial) 
one and a convected (material) coordinate system. The relations between the 
time derivative and the covariant derivatives (gradients) have been developed 
for curved space and motion, so that the velocity components supply the 
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connection between the equations of motion and the time rate of change of the 
metric and curvature tensors. 

The metric tensor (time rate of change) in the convected material coordi- 
nate system is linearly decomposed into elastic and plastic parts. In this 
formulation, a yield function is assumed, which is dependent on the rate of 
change of stress, metric, temperature, and a set of internal variables. 
Moreover, a hypoelastic law was chosen to describe the thermoelastic part of 
the deformation. 

A time and temperature dependent viscoplastic model was formulated in 
this convected material system to account for finite strains and rotations. 
The history and temperature dependence were incorporated through the introduc- 
tion of internal variables. The choice of these variables, as well as their 
evolution, was motivated by phenomenological thermodynamic considerations. 

The nonisothermal elastic-viscoplastic deformation process was described 
completely by "thermodynamic state" equations. Most investigators (in the 
area of viscoplasticity) employ plastic strains as state variables. Our study 
shows that, in general, use of plastic strains as state variables may lead to 
inconsistencies with regard to thermodynamic considerations. Furthermore, the 
approach and formulation employed by all previous investigators lead to the 
condition that all plastic work is completely dissipated. This, however, is 
in contradiction with experimental evidence, from which it emerges that part 
of the plastic work is used for producing residual stresses in the lattice, 
which, when phenomenologically considered, causes hardening. Both limitations 
are not present in our formulation, because of the inclusion of the "thermody- 
namic state" equations. 

The obtained complete rate field equations consist of the principles of 
the rate of the virtual power and the rate of conservation of energy, of the 
constitutive relations, and of boundary and initial conditions. These formu- 
lations provide a sound basis for the formulation of the adopted finite 
element solution procedures. 

The derived shell theory, in the least restricted form, before any 
simplifying assumptions are imposed, has the following desirable features: 

(a) The two-dimensional, impulse-integral form of the equations of 
motion and the Second Law of Thermodynamics (Clausius-Duhem inequal- 
ity) for a shell follow naturally and exactly from their three-di- 
mensional counterparts. 

(b) Unique and concrete definitions of shell variables such as stress 
resultants and couples, rate of deformation, spin and entropy 
resultants can be obtained in terms of weighted integrals of the 
three-dimensional quantities through the thickness. 

(c) There are no series expansions in the thickness direction. 

(d) There is no need for making use of the Kirchhoff Hypotheses in the 
kinematics. 

(e) All, approximations can be postponed until the descretization process 
of the integral forms of the First Law of Thermodynamics 

(f) A by-product of the descent from three-dimensional theory is that 

the two-dimensional temperature field (that emerges) is not a 
through-the-thickness average, but a true point by point 

distribution. This is contrary to what one finds in the literature 
concerning thermal stresses in the shell. 


To develop geometrically nonlinear, doubly curved finite shell elements 
the basic equations of nonlinear shell theories have to be transferred into 
the finite element model. As these equations in general are written in tensor 
notation, their implementation into the finite element matrix formulation 
requires considerable effort. 

The nonlinear element matrices are directly derived from the incremental- 
ly formulated nonlinear shell equations, by using a tensor-oriented procedure. 
The classical thin shell theory based on the Kirchoff-Love hypotheses (Formu- 
lation D in Appendix A) was employed for this purpose. For this formulation, 
we are using the "natural" degrees of freedom per mid-surface shell node: 
three incremental velocities and the rates of rotations about the material 
coordinates in a mixed form. 

A description of the developed element and related finite element code 
are given in Appendix B. This exposition provides information concerning the 
formulation, the finite element and how it is employed in the solution of 
shell-like configurations. Moreover a complete description including program 
flow chart, listing, input instructions to the user and explanation of output 
are also included in Appendix B. 

The quasi -linear nature of the principle of the rate of virtual power 
suggests the adoption of an incremental approach to numerical integration with 
respect to time. The availability of the field formulation provides assurance 
of the completeness of the incremental equations and allows the use of any 
convenient procedure for spatial integration over the domain V. In the 
present instance, the choice has been made in favor of a simple first order 
expansion in time for the construction of incremental solutions from the 
results of finite element spatial integration of the governing equations. 

The procedure employed permits the rates of the field formulation to be 
interpreted as increments in the numerical solution. This is particularly 
convenient for the construction of incremental boundary condition histories. 

Even under the condition of static external loads and slowly growing 
creep effects, the presence of snap-through buckling makes the inertial 
effects significant. In dynamic analyses, the applied body forces include 
inertial forces. Assuming that the mass of the body considered is preserved, 
the mass matrix can be evaluated prior to the time integration using the 
initial configuration. 

Finite element solution of any boundary-value problem involves the 
solution of the equilibrium equations (global) together with the constitutive 
equations (local). Both sets of equations are solved simultaneously in a step 
by step manner. The incremental form of the global and local equations can be 
achieved by taking the integration over the incremental time step t=t. +1 -t.. 
The rectangular rule has been applied to execute the resulting time integra- 
tion. 


Clearly, the numerical solution involves iteration. A simplified version 
of the Riks-Wempner constant-arc-length method has been utilized. This itera- 
tion procedure which is a generalization of the displacement control method 
also allows to trace the nonlinear response beyond bifurcation points. In 
contrast to the conventional Newton-Raphson techniques, the iteration of the 
method takes place in the velocity and load rate space. The load step of the 
first solution in each increment is limited by controlling the length ds of 
the tangent. Either the length is kept constant in each step or it is adapted 
to the characteristics of the solution. In each step the triangular-size 
stiffness matrix has to be checked for negative diagonal terms, indicating 
that a critical point is reached. 

One of the most challenging aspects of finite strain formulations is to 
locate analytical solution with which to compare the proposed formulation. 
Typically, as a first problem, a large strain uniaxial test case was analyzed. 
The case considered examines the rate-dependent plastic response of a bar to a 
deformation history that includes segments of loading, unloading, and reload- 
ing, each occurring at varying strain and temperature rates. Moreover, it was 
shown that the proposed formulation generates no strain energy under a pure 
rigid body rotation. These are surely important demonstrations but they only 
represent a partial test because the principal stretch directions remain 
constant. Finally, a problem which was discussed by Nagtegaal and de Jong, 
and others too, as a problem which demonstrates limitations of the constitu- 
tive models in many strain formulation,, is the Couette flow problem. This 
problem is solved as a third example. The results of these test problems show 
that: 

The formulation can accommodate very large strains and rotations. 

The formulation does not display the oscillatory behavior in the 
stresses of the Couette flow problem. 

The model incorporates the simplification associated with 
rate-insensitive elastic response without losing the ability to 
model rate temperature dependent yield strength and plasticity. 

The problem of buckling of shallow arches under transient 

thermomechanical load was investigated next. The analysis was performed with 
the aid of 24 paralinear isoparametric elements. The paral inear isoparametric 
element is such that the thickness is small compared to other dimensions. The 
characteristics of the element are defined by the geometry and interpolation 
functions, which are linear in the thickness direction and parabolic in the 
longitudinal direction. Consequently, the element allows for shear strain 

energy since normals to a mid-surface before deformation remain straight, but 
not necessarily normal to the midsurface after deformation. 

The developed solution scheme is capable of predicting response which 
includes pre- and post-buckling with thermal creep and plastic effects. The 
solution procedure was demonstrated through several examples which include 
both creep and snap-through behavior. 

The last set of problems which are under investigation consists of creep 
or thermal buckling, with plastic effects, of shells of revolution. 



In addition, following a more traditional approach, a method was devel- 
oped for bounding the response (solution) of bars and beams of (linear) 
visoelastic material behavior, based on nonlinear kinematic relations. 

In connection with the progress to date, two papers were published by the 
AIAA Journal in 1986 and 1987. Moreover, a paper entitled, "Non-Isothermal 
Elastoviscoplastic Analysis of Planar Curved Beams" was presented at the 3rd 
Symposium on Nonlinear Constitutive Relations for High-Temperature Applica- 
tions, held at the University of Akron, on June 11-13, 1986. A descriptive 
abstract of this paper was published in the meeting proceedings and the full 
paper appeared in a special publication (NASA CP 10010). Copies of the above 
have been sent to the sponsor. 

In addition, the two papers presented at the 28th AIAA/ASME/ASCE/AHS SDM 
Conference and published in the proceedings of this conference, have been 
accepted for publication; both by the AIAA Journal. These two papers deal 
with applications to snap-through and creep buckling of bars and arches. Most 
of this work was also presented at the NASA-Lewis Conference on Structural 
Integrity and Durability of Reusable Space Propulsion Systems on May 1987 in 
Cleveland, Copies of these papers will be forwarded to the Sponsor, as soon as 
they appear in print. 

In connection with the more traditional approach a paper accepted for 
presentation and for publication in the Proceedings of the special Symposium 
on Constitutive Equations at the ASME Winter Annual Meeting, Chicago, IL., 
November 28 - December 2, 1988. The title of the paper is "Creep Analysis of 
Beams and Arches Based on a Hereditary Visco- Elastic- Plastic Constitutive 
Law". Copies of this paper have been forwarded to the sponsor in December, 
1988 .. 


Moreover, one paper was published in the the AIAA/ASME/. . ./AHS 30th SDM 
Conference. A copy of this paper is attached, herewith. 

Finally, a finite element has been developed and it is currently being 
tested for a less restricted shell formulation (Formulation C; see Appendix 
A). A description of this newly developed element and the related finite 
element code will be submitted to the sponsor as soon as the testing is 
completed and reliable results have been obtained. 


FUTURE TASKS 


The main thrust of the additional tasks is to develop a finite element 
and select a code, which will be made available to all users and which will be 
based on the most general (but practical) nonlinear shell formulation possible 
and nonlinear constitutive relations to predict the response of shell-like 
structures, when subjected to time-dependent thermomechanical loads with large 
excursions. 

This should be completed in the next six months and the complete package 
will be delivered to the sponsor. 


APPENDIX A 


The various shell theory approximations (formulations) are based on the 
use of certain simplifying assumptions regarding the geometry and kinematics 
of the shell configuration. 

These are: 

Assumption I : The material points which are on the normal to the refer- 

ence surface before deformation will be on the same normal after deforma- 
tion. 

Assumption II : The shell is sufficiently thin so that we can assume 

linear dependence of the position of any material point (in the deformed 
state) to the normal (to the reference surface) coordinate (in the 
deformed state). The linear dependence can easily be changed to parabol- 
ic, cubic, or any desired degree of approximation. 

Assumption III : The rate of change of the velocity gradients with 

respect to in-plan coordinates on the two boundary shell surfaces is 
negligibly small . 

Assumption IV : The rate of change of the distance of a material from the 

reference surface is negligibly small. 

On the basis of the above four simplifying assumptions, several formula- 
tions result, for the analysis of thin shells. 

This formulation makes use of Assumption I, only. 

This formulation employs Assumptions I and II. 

This formulation employs Assumption I, II and III. 

This is the classical thin shell theory based on the 
Kirchhoff-Love hypotheses of Assumptions I, II, III, IV, 
as applied to large deformation theory. 

These formulations are arranged in such a manner that we move from the least 
restrictive (A) to the most restrictive (D). 

In addition to this a fifth formulation (E) can easily be devised and 
this formulation in terms of order of restriction is similar to Formulation A. 
Formulation E makes use of Assumption II only. 

APPENDIX B FINITE ELEMENT AND RELATED CODE FOR FORMULATION D 

In this Appendix, a description of the developed finite element and the 
related finite element code is presented. First, the essentials of the 
element are described and then the complete solution procedure is presented 
with sufficient detail. This includes a flow chart, a line by line listing of 
the computer program, input data information, and explanation of the output. 


Formulation A 
Formulation B 
Formulation C 
Formulation D 
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6,1. THE SHELL ELEMENT 

A brief description highlighting the essential features of 
the shell element development and the related code used in this 
work is given here. 

In order to derive discrete algorithm based on the finite 
element displacement method we approximate the velocity field by 
index-oriented notation, which allows the separate representation 
of the shape functions (the specific expression depends on the 
decided upon degress - of - freedom, Lagrangian, Hermitian, etc.) 
for the tangential velocities and for the normal velocity 0* . 

VJi= V/\ - V" ^ (1) 

VJj Vav * V* (2) 

Upper indices imply the columns, lower indices the rows of a 
matrix expression, and the summation is carried out spanning the 
number of degrees-of-freedom. V M and V M represent, therefore, the 
vector of nodal velocities by the row and by the column respec- 
tively. We get the shape functions for the partial derivatives 
of the velocity shape functions } %t* : 


a< lP =(vr2!) (P V/»=- Vm 

'Vr(' r »),p v *' ^ 


The main idea of this formulation is the development of 
shape functions for further mechanical and thermal variables by 
the application of well-known tensor procedure on the basic shape 
functions (3) and (4). Taking, for example, the covariant 
derivative 

Ik; £ a ( 5 ) 

and inserting (3) and (1) , we can define the shape function of 
the covariant' derivative: 

vJV*j e> * ( > ~ P oip. ) Vm - VJV& Vm (6) 

M 

In the same way we receive the shape function of the 
internal variables, for example the rate of deformation and the 
spin tensors: 

i ( r: it * »£* - ^ )Va = a; p) Vv <7) 

v> y — J 




(8) 
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The same procedure is now applied to the shift tensor 


+ st < 9 > 

which is responsible for the "exact" distributions over the 
thickness. 

Finally the shape functions of the internal forces and 
temperature variables can be derived from the shape functions of 
the rate of deformation and spin tensor via the constitutive 
relations; for example: 

/•“« = p u-,%, V- ■ M V* (io, 

All these expression are now substituted into the rate of 
the first law of thermodynamics to obtain the element "stiffness" 
equations. The developed element matrices are implemented into 
the global relation of the complete shell structure by standard 
assemblage process considering incidence and boundary conditions. 

The present curvilinear formulation of the element enables 
the precise description of geometry, external loads and tempera- 
tures and the fulfillment of the convergence criteria, while the 
rigid body motion condition can only be satisfied in an approxi- 
mate manner. The tensor oriented formulation renders the 
optional use of various shape functions for the tangential 
velocities v and the normal velocity V 3 . 

The shell element which have been used up to today is based 
on the bicubic Hermite polynomial with 4x12 generalized velocit- 
ies and 4 temperatures. Numerical integration spanning the 
element domain was applied (16 points of integration) , whereby 
area and boundary integrals were replaced by double integration 
with respect to the curvilinear Q 4- - coordinates. 
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